Estimation ========== RISE estimates models by three approaches through a single entry point: .. list-table:: :header-rows: 1 :widths: 22 32 46 * - Approach - When to use it - How to call * - **Bayesian** - Informative priors, fragile identification, or model comparison via posterior model probabilities. - ``estimate(m, estim_priors = priors)`` * - **Maximum likelihood** - Data-only inference, comparison with Bayesian, or classical inference on a non-DSGE shape. - ``estimate(m, estim_priors = boundsOnly)`` (uniform-bounds priors) * - **Indirect inference** - Likelihood intractable; match specific moments, IRFs, or auxiliary statistics. - ``indirect_inference(m, myobjective, estim_priors = priors)`` All three share the same mode-finding optimisers, the same restriction machinery, the same Hessian / standard-error pipeline, and the same visualisation tooling. The differences are entirely in the objective and in which post-estimation outputs are meaningful. This chapter is the modern overview of the estimation surface. The in-depth treatment is split into the sub-chapters below (restrictions, priors, posterior maximization, posterior simulation, marginal data density, processing posterior draws, indirect inference), mirroring the legacy structure. .. toctree:: :maxdepth: 2 Estimation restrictions Prior distributions and implementation Posterior maximization Posterior simulation Indirect Inference .. contents:: :local: :depth: 2 The estimation pipeline at a glance ------------------------------------ A complete estimation run is four steps:: % 1. Bind data and priors to a parsed model m = set(m, data = mydata, estim_priors = priors); % 2. Find the posterior mode m_est = estimate(m); % 3. (Optional) sample the posterior [target, x0, lb, ub, SIG] = pull_objective(m_est); obj = rsamplers.rwmh(target, x0, lb, ub, ... struct(N = 2000, tunedCov = SIG)); results = sample(obj); % 4. (Optional) diagnose the chain mc = mcmc(results); mc.summary(); mc.traceplot('phi_pi'); Dispatch on properties, not class ---------------------------------- A modern-toolbox distinctive worth flagging up front: the estimator chooses its likelihood path based on **properties of the data and the model**, never the model's class. Concretely: * If the data has missing observations or unobserved states, the estimator runs a Kalman filter (or the regime-conditional Kalman filter, if ``h > 1``). * If the data are complete and the model admits a closed-form posterior moment (notably BVAR with a conjugate prior), the estimator short-circuits to the closed form. * For nonlinear models the estimator runs the appropriate nonlinear filter (regime-switching unscented or particle filter). User code that constructs a model from a new shape gets correct estimation behavior for free -- there is no estimator branch keyed on class name. See :doc:`../Architecture/Modern architecture` for the broader pattern. Priors ------- A prior is a struct whose fields are parameter names and whose values are cell arrays describing the prior. RISE accepts **five parametrisations** for a single-parameter prior, plus a Dirichlet form for transition probabilities and a rich *endogenous priors* facility for behavioral constraints. Common abbreviations used below: * ``start`` -- initial value for the optimizer. * ``mean`` -- mean of the prior distribution. * ``sd`` -- standard deviation. * ``lb``, ``ub`` -- absolute bounds (hard truncation). * ``lq``, ``uq`` -- lower and upper quantiles. * ``hpp1``, ``hpp2``, ``hpp3``, ``hpp4`` -- raw hyperparameters. * ``p`` -- the struct of priors; ``alpha`` -- a parameter name. Uniform priors ~~~~~~~~~~~~~~~ :: p.alpha = {start, lb, ub}; Internally rewritten to ``{start, lb, ub, 'uniform(1)'}``. If you set every estimated parameter this way, you have done **maximum likelihood** (see Ireland, 2004) -- the posterior reduces to the likelihood up to a constant. Means and standard deviations ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ :: p.alpha = {start, mean, sd, 'distribution'}; With hard truncation:: p.alpha = {start, mean, sd, 'distribution', lb, ub}; With more hyperparameters (e.g. generalized beta, t-student, truncated normal):: p.alpha = {start, mean, sd, hpp3, hpp4, 'distribution'}; p.alpha = {start, mean, sd, hpp3, hpp4, 'distribution', lb, ub}; Quantiles ~~~~~~~~~~ Means and standard deviations do not always exist (e.g. inverse gamma with certain hyperparameters; Cauchy). Quantile-based parametrisation is the universal workaround:: p.alpha = {start, lq, uq, 'distribution(probability)'}; The ``(probability)`` suffix is the mass that should fall between ``lq`` and ``uq`` -- typically ``0.9``. With truncation and additional hyperparameters:: p.alpha = {start, lq, uq, 'distribution(probability)', lb, ub}; p.alpha = {start, lq, uq, hpp3, hpp4, 'distribution(probability)'}; p.alpha = {start, lq, uq, hpp3, hpp4, 'distribution(probability)', lb, ub}; Direct hyperparameters ~~~~~~~~~~~~~~~~~~~~~~~ The previous parametrisations are convenience wrappers around the distribution's hyperparameters. You may set them directly:: p.alpha = {start, hpp1, hpp2, 'distribution(hpp)'}; with truncation and extra hyperparameters as for the other forms. When using this parametrisation, plot the distribution and inspect it carefully before estimating; it is the most error-prone way. User-defined priors via PDF ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ If you need a prior RISE does not support, supply the **PDF** (not the log-PDF) as a function handle plus hard bounds:: p.alpha = {start, lb, ub, function_handle}; RISE constructs the CDF, the inverse CDF, and everything else internally. The t-student special case ~~~~~~~~~~~~~~~~~~~~~~~~~~~ The t-student's third hyperparameter is the degrees of freedom; there is no fourth, so leave it as ``nan`` or empty:: p.alpha = {start, 3, 7, 10, nan, 't(.9)'}; Dirichlet priors for transition matrices ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For a ``n x n`` transition matrix whose entries are named ``pname(i,j)``, declare one Dirichlet per row:: p.dirichlet_1 = {stdev(1,1), pname(1,2), mean(1,2), ... pname(1,3), mean(1,3), ..., pname(1,n), mean(1,n)}; with ``mean(1,1) + mean(1,2) + ... + mean(1,n) = 1`` (the diagonal mean is implied) and ``stdev(1,1)`` the standard deviation of the diagonal element. The ``_1``, ``_2``, ..., ``_n`` suffix is just a label -- the order is irrelevant and you only need one Dirichlet per row you actually want to put a prior on. Endogenous (system / property / behavioral) priors ---------------------------------------------------- Rather than placing priors on parameters, you may place priors on **model-implied behaviors**: covariances, IRFs, variance and historical decompositions, first-order solution coefficients, regime probabilities at specific dates, aggregates of any of the above. The literature calls these *endogenous priors* (Del Negro & Schorfheide, 2008; Christiano et al., 2011), *system priors* (Andrle & Benes, 2013; Andrle & Plasil, 2017, 2018), *property priors* (Tetlow), or *behavioral priors*. The general syntax mirrors the parameter-prior cell shapes, with the parameter name replaced by an **expression string**:: {'expression', mean, std, 'distribution'} {'expression', lb, ub} {'expression', lq, uq, 'distribution(probability)'} {'expression', hpp1, hpp2, 'distribution(hpp)'} The expression syntax covers seven classes; regime arguments default to ``1`` in switching models when omitted. Covariance and correlation ~~~~~~~~~~~~~~~~~~~~~~~~~~~ :: 'cov{var1, var2, lag}' 'cov{var1, var2, lag, regime}' 'corr{var1, var2, lag}' 'corr{var1, var2, lag, regime}' ``'cov{C, Y, 1, 2}'`` is the covariance between ``C`` and ``Y`` at lag 1 in regime 2. Impulse responses ~~~~~~~~~~~~~~~~~~ :: 'irf{variable, shock, horizon}' 'irf{variable, shock, horizon, regime}' ``'irf{C, EA, 3, 2}'`` is the response of ``C`` to shock ``EA`` at horizon 3 in regime 2. Variance decomposition ~~~~~~~~~~~~~~~~~~~~~~~ Finite horizon:: 'vd{variable, shock, horizon}' 'vd{variable, shock, horizon, regime}' Long run (set ``horizon = inf``):: 'vd{variable, shock, inf}' 'vd{variable, shock, inf, regime}' Historical decomposition ~~~~~~~~~~~~~~~~~~~~~~~~~ :: 'hd{variable, shock, horizon}' First-order solution coefficients ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Impact of a state variable on an endogenous variable:: 'Tz{variable, stateVariable}' 'Tz{variable, stateVariable, regime}' Examples: * ``'Tz{C, C{-1}}'`` -- impact of ``C{-1}`` on ``C`` in regime 1. * ``'Tz{C, @sig}'`` -- impact of the perturbation parameter. * ``'Tz{C, EA}'`` -- impact of shock ``EA``. * ``'Tz{C, EA{+3}, 2}'`` -- expected impact of ``EA`` three periods ahead, in regime 2. Filtered / updated / smoothed probabilities ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Regime probabilities:: 'mvar{f, regime, date_range}' % filtered 'mvar{u, regime, date_range}' % updated 'mvar{s, regime, date_range}' % smoothed State combinations across chains:: 'mvar{f, state1 & state2 & ... & stateN, date_range}' ``'mvar{f, vol_2 & pol_1, rq(2020,3)}'`` -- filtered probability of being simultaneously in ``vol_2`` and ``pol_1`` in 2020Q3. Aggregations and discontinuous ranges ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ``sum``, ``max``, ``min`` over a date or horizon range:: 'sum(irf{C, EA, 1:40, 3})' 'max(mvar{f, regime_2, rq(1970,3):rq(2024,3)})' Discontinuous ranges:: 'mvar{f, regime_2, [rq(1970,3):rq(1975,3), rq(1978,3), rq(1980,3):rq(1985,4)]}' User-defined endogenous priors ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ A custom function on the model:: {@(model) my_function(model), mean, std, 'distribution'} {@(model) my_function(model), lb, ub} For advanced cases where the prior needs filtration information, the function handle, called with **one** input (the model object), returns a struct:: v = struct('priors', prior_cell, ... 'kf_filtering_level', 0); with ``kf_filtering_level`` in: * ``0`` -- no filtering information required; * ``1`` -- one-step-ahead (filtered); * ``2`` -- updated; * ``3`` -- smoothed. When called with **two** inputs (model + filtering struct), the function returns a vector of numbers evaluating the priors. List of supported prior distributions -------------------------------------- The complete catalogue, with their supports. Names with ``_pdf`` suffix (e.g. ``beta_pdf``) are accepted as synonyms. .. list-table:: :header-rows: 1 :widths: 25 25 50 * - Name - Support - Notes * - ``arcsine`` - :math:`[0, 1]` - * - ``beta`` - :math:`[0, 1]` - Conjugate for Bernoulli/binomial. * - ``dirichlet`` - simplex - For transition-matrix rows; see above. * - ``betaprime`` - :math:`[0, \infty)` - * - ``cauchy`` - :math:`(-\infty, \infty)` - No mean; quantile parametrisation only. * - ``exponential`` - :math:`[0, \infty)` - * - ``folded_normal`` - :math:`[0, \infty)` - * - ``gamma`` - :math:`[0, \infty)` - * - ``gammatau`` - :math:`[0, \infty)` - Generalised gamma-tau. * - ``gumbel`` - :math:`(-\infty, \infty)` - * - ``half_normal`` - :math:`[0, \infty)` - * - ``igamma`` (= ``igamma2``) - :math:`[0, \infty)` - Use for standard deviations. * - ``igamma1`` - :math:`[0, \infty)` - Inverse-gamma, shape/scale parametrisation. Use for standard deviations. * - ``igamma2`` - :math:`[0, \infty)` - Use for standard deviations. * - ``inverse_gaussian`` - :math:`(0, \infty)` - Use for standard deviations. * - ``kumaraswamy`` - :math:`(0, 1)` - Alternative to ``beta`` on the unit interval. * - ``laplace`` - :math:`(-\infty, \infty)` - * - ``left_triang`` - :math:`[a, b]` - Skewed to the left. * - ``levy`` - :math:`[\mu, \infty)` - * - ``logistic`` - :math:`(-\infty, \infty)` - * - ``lognormal`` - :math:`(0, \infty)` - * - ``normal`` - :math:`(-\infty, \infty)` - * - ``pareto`` - :math:`[x_m, \infty)` with :math:`x_m > 0` - * - ``raised_cosine`` - :math:`[\mu - s, \mu + s]` with :math:`s > 0` - * - ``sichisq`` - :math:`(0, \infty)` - Scaled inverse :math:`\chi^2`. Use for standard deviations. * - ``t`` - :math:`(-\infty, \infty)` - DOF is the third hyperparameter; no fourth (set ``nan``). * - ``truncated_normal`` - :math:`[\underline{\theta}, \overline{\theta}]` - * - ``uniform`` - :math:`[\underline{\theta}, \overline{\theta}]` - Implied when only the three-slot form is used. * - ``weibull`` - :math:`[0, \infty)` - Visualizing priors (and posteriors) ------------------------------------ ``rdist.plot`` plots the priors and, when given an ``@mcmc`` object, overlays the posterior densities:: rdist.plot(priors); rdist.plot(priors, plotOpts); rdist.plot(priors, plotOpts, mc); rdist.plot(priors, plotOpts, mc, 'linewidth', 2); ``plotOpts`` is a struct with the following fields: .. list-table:: :header-rows: 1 :widths: 25 75 * - Field - Effect * - ``r0c0`` - Vector ``[nrows, ncols]`` of the figure grid. Default ``[3, 3]``. * - ``prior_trunc`` - Density-tail truncation. Default ``1e-3``. * - ``scale_IG_trunc`` - Correction of ``prior_trunc`` for inverse-gamma families. Default ``25``. * - ``npoints`` - Number of points for the density evaluation. Default ``1000``. * - ``do`` - Struct with logical fields ``priors``, ``posteriors`` -- whether to draw each. * - ``par_list`` - Cell array of parameter names to restrict the plot to. Example:: priors = struct(); priors.beta = {.5, 95, 5, 'beta(hpp)'}; priors.delta = {.5, 40, 860, 'beta(hpp)'}; priors.alpha = {.5, 360, 590, 'gamma(hpp)'}; priors.rhoz = {.5, 80, 20, 'beta(hpp)'}; priors.sdez = {0.1, 0.01, 0.01, 'igamma(hpp)'}; plotOpts = struct(); plotOpts.prior_trunc = 2.6e-3; rdist.plot(priors, plotOpts, [], struct('linewidth', 2)); Where priors live in the toolbox --------------------------------- By convention, calibration and priors live together in a companion ``mymodel_est_params.m`` file with the signature ``function [p, priors] = mymodel_est_params()``. See :doc:`../ModelShapes/DSGE/Model file language` for why these stay out of the ``.rs`` file. Switching parameters use their flat names with chain and regime appended (``phi_pi_mon_1``, ``mon_tp_1_2``, ...). Transition probabilities under a Dirichlet prior use the ``dirichlet_*`` naming pattern shown above. Maximum likelihood (no priors) ------------------------------- To estimate by maximum likelihood, supply only **bounds** for each parameter -- no mean, no distribution:: priors.alpha = {start, lb, ub}; priors.rho = {start, lb, ub}; priors.sigma_a = {start, lb, ub}; m_est = estimate(m, data = mydata, estim_priors = priors); Internally RISE rewrites each entry as ``{start, lb, ub, 'uniform(1)'}``; the prior contribution becomes a constant on the box and the posterior reduces to the (truncated) likelihood. :cite:`Ireland2004` is the canonical DSGE reference for this practice. What "no prior" means in practice: * The mode is the **maximum-likelihood estimator**, not a posterior mode. * The Hessian-derived covariance is the standard ML covariance :math:`H = -\partial^2 \log L / \partial \theta^2` at the mode. * Bayesian-only outputs are reinterpreted: the marginal data density reverts to the integrated likelihood over the box, not the Bayesian posterior probability of the model. * Posterior simulation still runs -- the chain explores the likelihood-times-indicator (uniform prior). Use it for profile-likelihood-style sensitivity rather than model comparison. Mixing approaches per parameter is allowed: some parameters can carry informative priors and others can be left with bounds-only (uniform). The dimension of the prior contribution adapts automatically. When the likelihood is intractable -- nonlinear filters that fail on your sample, or quantities of interest that depend on simulated paths -- consider :doc:`indirect inference ` instead. The optimiser suite and restriction machinery are the same. Mode finding ------------- With data and priors bound, ``estimate`` runs an optimizer on the log-posterior:: m_est = estimate(m); The default optimizer is ``fmincon``. Pick a different one with the ``optimizer`` option:: m_est = estimate(m, optimizer = 'wcsminwel'); To pass options through to the optimizer, use the cell form ``{name, opts}``:: opt = optimset('fmincon'); opt.MaxIter = 2000; opt.TolFun = 1e-8; opt.Display = 'iter'; m_est = estimate(m, optimizer = {'fmincon', opt}); Bundled optimizers ~~~~~~~~~~~~~~~~~~~ .. list-table:: :header-rows: 1 :widths: 18 22 14 12 34 * - Optimizer - Family - Scope - Parallel - Typical use * - ``fmincon`` - SQP / interior-point - local - via ``UseParallel`` - default; smooth likelihoods. * - ``fminunc_bnd`` - quasi-Newton (BFGS) - local - no - smooth, only-bound constraints. * - ``fminsearch_bnd`` - Nelder-Mead - local - no - very small problems; non-smooth. * - ``patternsearch`` - direct search (GPS / MADS) - local - via Global Opt Toolbox - non-smooth / noisy objectives. * - ``wcsminwel`` - Sims's csminwel + restart - local - no - robust on awkward DSGE surfaces. * - ``wnewrat`` - Newton with numerical Hessian - local - no - mode polishing after a global pass. * - ``wgmhmaxlik`` - MH-driven climb - local + exploration - no - pre-conditioner for multi-modal surfaces. * - ``bee_gate`` - artificial bee colony - global (population) - yes - older global default. * - ``rise_lshade``, ``rise_jso``, ``rise_cma_es``, ``rise_de``, ``rise_abc``, ``rise_aco``, ``rise_bbo``, ``rise_agsk``, ``rise_mads`` - ``+globalopt`` family - global - yes (per algorithm) - the modern global search menu. * - ``sa.gmhmaxlik`` - simulated annealing - global - no - robust, slow; legacy. * - ``blockwise_optimization`` - dispatcher - same as inner - same as inner - wraps any optimizer to optimize block-by-block; triggered automatically when ``estim_blocks`` is non-empty. Common options accepted across optimizers ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. list-table:: :header-rows: 1 :widths: 25 60 15 * - Field - Meaning - Default * - ``MaxIter`` - iteration cap - ``1000`` * - ``MaxFunEvals`` - total function-eval cap - ``inf`` * - ``MaxTime`` - wall-clock cap (seconds) - ``inf`` * - ``TolFun`` - objective tolerance - ``1e-6`` * - ``TolX`` - step tolerance - ``1e-6`` * - ``Display`` - ``'iter'`` / ``'on'`` / ``'off'`` - ``'off'`` * - ``MaxNodes`` - population size (population methods) - algorithm-specific * - ``ObjectiveLimit`` - early-stop floor - ``-inf`` Algorithm-specific hyperparameters (``F``, ``CR`` for ``rise_de``; ``sigma0`` for ``rise_cma_es``; ``limit`` for ``rise_abc``; ...) are documented in legacy *Stochastic Global Optimization*. Running several optimizers in turn ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ A common strategy is global search to locate the basin, then a fast local method to polish:: m_est = estimate(m, ... optimizer = {'rise_lshade', 'wcsminwel'}); Each optimizer starts from the previous one's best point. User-defined optimizers ~~~~~~~~~~~~~~~~~~~~~~~~ ``optimizer`` may also be a function name or a function handle implementing:: [xfinal, ffinal, exitflag, H] = ... myOptimizer(fh, x0, lb, ub, options, varargin) ``fh`` is the negative log-posterior to **minimize**; RISE guarantees ``lb <= x0 <= ub`` and a clean options struct. Outputs: ``xfinal`` (column vector in bounds), ``ffinal`` (scalar value ``fh(xfinal)``), ``exitflag`` (MATLAB convention: ``1`` = converged, ``0`` = budget exhausted, negative = failed), ``H`` (``n x n`` Hessian at ``xfinal``, or ``[]`` if not available -- the estimator falls back to a finite-difference pass). Manual stopping is honoured by every ``+globalopt`` optimizer: drop a file ``ManualStoppingFile.txt`` in the working directory and the next iteration boundary returns the current best. What ``estimate`` puts on the model ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ After ``estimate`` returns:: m_est.estimation.posterior_maximization.mode m_est.estimation.posterior_maximization.mode_stdev m_est.estimation.posterior_maximization.hessian m_est.estimation.posterior_maximization.vcov m_est.estimation.posterior_maximization.log_post m_est.estimation.posterior_maximization.log_lik m_est.estimation.posterior_maximization.log_prior m_est.estimation.posterior_maximization.log_marginal_data_density_laplace Parameter names corresponding to the mode vector are ``fieldnames(m_est.estimation.priors)`` in declaration order. Posterior simulation --------------------- The mode is rarely the end. To explore the posterior, run an MCMC sampler. The modern API is **class-based** -- each sampler is a class in ``+rsamplers`` constructed from a target function and bounds, then sampled with the ``sample`` method. ``pull_objective`` lifts a sampler-ready objective out of the estimated model:: [target, x0, lb, ub, SIG] = pull_objective(m_est); * ``target`` is the log-posterior at a parameter vector. It is the value to **maximize** (note the sign: this differs from the legacy ``sampler_rwmh`` shortcut, which expected a value to *minimize*). * ``x0``, ``lb``, ``ub`` are the starting point and bounds. * ``SIG`` is the Hessian-derived covariance, a sensible proposal covariance. Available samplers ~~~~~~~~~~~~~~~~~~~ .. list-table:: :header-rows: 1 :widths: 28 72 * - Class - Algorithm * - ``rsamplers.rwmh`` - Random-walk Metropolis-Hastings. Defaults: ``c=1``, ``proposal='normal'``, ``tunedCov`` from the mode Hessian. * - ``rsamplers.imh`` - Independent Metropolis-Hastings. * - ``rsamplers.apt`` - Adaptive parallel tempering. Useful for multi-modal posteriors. * - ``rsamplers.slice`` - Slice sampler. No proposal covariance to tune. * - ``rsamplers.usrsmplr`` - User-defined sampler hook. Options shared by every sampler (inherited from ``rsamplers.rsampler``): .. list-table:: :header-rows: 1 :widths: 25 15 60 * - Field - Default - Role * - ``N`` - ``2000`` - Draws per chain. * - ``nchain`` - ``1`` - Number of chains. * - ``burnin`` - ``0`` - Burn-in draws (discarded). * - ``thinning`` - ``1`` - Keep every ``thinning``-th draw. * - ``MaxTime`` - ``inf`` - Wall-clock cap (seconds). * - ``MaxFunEvals`` - ``inf`` - Total function-eval cap. Constructing and running ~~~~~~~~~~~~~~~~~~~~~~~~~ :: opts = struct( ... N = 2000, ... nchain = 2, ... burnin = 500, ... c = 0.5, ... tunedCov = SIG); obj = rsamplers.rwmh(target, x0, lb, ub, opts); results = sample(obj); ``results`` is a cell of length ``nchain``. Each ``results{k}`` carries ``pop`` (struct array of draws with ``.x`` (parameter vector) and ``.f`` (log-posterior), ``stats``, ``SIG``, ``c``, ``x0``, ``best``. Extract draws as a matrix:: draws = [results{1}.pop.x]; % nparams x N .. note:: The legacy toolbox exposed convenience shortcuts ``sampler_rwmh``, ``sampler_imh``, ``sampler_apt``, ``sampler_slice`` that expected a value to *minimize* and ran the sampler in one call. The modern path uses the class directly. The ``sampler_*`` shortcuts still work but should be migrated. User-defined sampler via ``rsamplers.usrsmplr`` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For samplers not in the bundle (e.g. DRAM, DA-DRAM, NUTS), write a wrapper that follows the RISE output convention: a cell of ``nchain`` structs, each with the ``pop``/``stats``/``SIG``/``c`` fields. The legacy *Posterior simulation* chapter contains a worked DRAM wrapper that illustrates the pattern; it remains canonical. Marginal data density ---------------------- For Bayesian model comparison, the ``mdd`` class wraps several estimators of :math:`\log p(y)`: * **Laplace** approximation -- automatic; available without a posterior simulation. * **Modified harmonic mean** (Geweke). * **Bridge sampling** (Meng-Wong). * **Mueller** (importance / reciprocal-importance). * **Chib-Jeliazkov** (auxiliary Gibbs). The Laplace approximation is cached on the estimated model:: log_mdd_laplace = ... m_est.estimation.posterior_maximization.log_marginal_data_density_laplace; For higher-quality estimators:: md = mdd(results, m_est); % uses the smc chain results log_mdd_mhm = md.modified_harmonic_mean(); log_mdd_bs = md.bridge_sampling(); log_mdd_cj = md.chib_jeliazkov(); Each method has its own option struct; see the legacy chapter for the per-estimator option surface. Chain diagnostics: ``@mcmc`` ---------------------------- The ``@mcmc`` class processes the output of ``sample`` and provides the standard convergence and quality diagnostics. The class itself has not changed between legacy and modern. Constructor:: mc = mcmc(results); mc = mcmc(results, chain_names); Properties .. list-table:: :header-rows: 1 :widths: 25 75 * - Property - Content * - ``pnames`` - Names of the estimated parameters. * - ``nchains`` - Number of chains. * - ``chain_names`` - Names assigned to the chains. * - ``npop`` - Number of population (parameter) draws. * - ``nparams`` - Number of parameters. * - ``draws`` - MCMC parameter draws stored as a 3D array (``nparams x npop x nchains``). * - ``psrf`` - Potential scale reduction factor (Gelman-Rubin). * - ``best`` - Best parameter draw. * - ``log_post`` - Log-posterior values for each draw. Diagnostic methods .. list-table:: :header-rows: 1 :widths: 30 70 * - Method - Effect * - ``traceplot(pname)`` - Per-chain trace plot of a parameter. * - ``densplot(pname)`` - Posterior density of a parameter. * - ``autocorrplot(pname)`` - Autocorrelation of a parameter's chain. * - ``meanplot(pname)`` - Running mean of a parameter's chain. * - ``scatterplot(pname1, pname2)`` - Joint scatter of two parameters. * - ``brooks_gelman()`` - Brooks-Gelman PSRF (requires ``nchain > 1``). * - ``plot_brooks_gelman()`` - Plot the PSRF. * - ``gelman_rubin()`` - Gelman-Rubin diagnostic. * - ``geweke(pname)`` - Geweke convergence diagnostic. * - ``raftery_lewis(pname)`` - Raftery-Lewis diagnostic. * - ``inefficiency_factors(pname)`` - Inefficiency factors. * - ``summary()`` - Compact tabular summary. * - ``posterior_parameter_statistics()`` - Posterior moments + credible intervals. .. warning:: ``mcmc`` is a **diagnostics class** constructed from existing draws. ``mcmc(m_est, ...)`` is not a sampler call and will error with "Brace indexing is not supported". To sample, use a ``rsamplers.*`` class; to diagnose, pass the result to ``mcmc``. Estimation restrictions ------------------------ Linear and nonlinear restrictions on the estimated parameters are declared at the model level via the ``estim_linear_restrictions`` and ``estim_nonlinear_restrictions`` options on ``set`` / ``estimate``. The general form is a cellstr of restrictions, one per row:: m = set(m, ... estim_linear_restrictions = { ... 'phi_pi - phi_y' ; ... 'mon_tp_1_2 - mon_tp_2_1' }); Each restriction is an expression that should equal zero. For inequalities, write them as differences with the appropriate sign convention. See the legacy chapter *Estimation restrictions* for the nonlinear and block-restriction extensions. The pipeline, skeletally ------------------------- :: % data, parameters, priors m = set(m, data = db, parameters = p, estim_priors = priors); % mode m_est = estimate(m); % a posterior sample [target, x0, lb, ub, SIG] = pull_objective(m_est); obj = rsamplers.rwmh(target, x0, lb, ub, ... struct(N = 2000, tunedCov = SIG)); results = sample(obj); % diagnose, compute the MDD mc = mcmc(results); log_mdd_mhm = mdd(results, m_est).modified_harmonic_mean();